{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a8b892c8",
   "metadata": {},
   "source": [
    "Printed and electronic copies of *Modeling and Simulation in Python* are available from [No Starch Press](https://nostarch.com/modeling-and-simulation-python) and [Bookshop.org](https://bookshop.org/p/books/modeling-and-simulation-in-python-allen-b-downey/17836697?ean=9781718502161) and [Amazon](https://amzn.to/3y9UxNb)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "collaborative-people",
   "metadata": {},
   "source": [
    "# Baseball"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "imported-table",
   "metadata": {
    "tags": []
   },
   "source": [
    "*Modeling and Simulation in Python*\n",
    "\n",
    "Copyright 2021 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-nc-sa/4.0/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "electoral-turkey",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# install Pint if necessary\n",
    "\n",
    "try:\n",
    "    from pint import UnitRegistry\n",
    "except ImportError:\n",
    "    !pip install pint\n",
    "    \n",
    "# import units\n",
    "from pint import UnitRegistry\n",
    "units = UnitRegistry()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "formal-context",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# download modsim.py if necessary\n",
    "\n",
    "from os.path import basename, exists\n",
    "\n",
    "def download(url):\n",
    "    filename = basename(url)\n",
    "    if not exists(filename):\n",
    "        from urllib.request import urlretrieve\n",
    "        local, _ = urlretrieve(url, filename)\n",
    "        print('Downloaded ' + local)\n",
    "    \n",
    "download('https://github.com/AllenDowney/ModSimPy/raw/master/' +\n",
    "         'modsim.py')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "progressive-typing",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# import functions from modsim\n",
    "\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "plastic-trigger",
   "metadata": {
    "tags": []
   },
   "source": [
    "This chapter is available as a Jupyter notebook where you can read the text, run the code, and work on the exercises. \n",
    "Click here to access the notebooks: <https://allendowney.github.io/ModSimPy/>."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "advanced-guidance",
   "metadata": {},
   "source": [
    "In the previous chapter we modeled an object moving in one dimension, with and without drag. Now let's move on to two dimensions, and baseball!\n",
    "\n",
    "In this chapter we model the flight of a baseball including the effect\n",
    "of air resistance. In the next chapter we use this model to solve an\n",
    "optimization problem."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "offensive-chinese",
   "metadata": {},
   "source": [
    "## Baseball\n",
    "\n",
    "To model the flight of a baseball, we have to make some decisions. To get started, we'll ignore any spin that might be on the ball, and the resulting Magnus force (see <http://modsimpy.com/magnus>). Under this assumption, the ball travels in a vertical plane, so we'll run simulations in two dimensions, rather than three.\n",
    "\n",
    "To model air resistance, we'll need the mass, frontal area, and drag\n",
    "coefficient of a baseball. Mass and diameter are easy to find (see\n",
    "<http://modsimpy.com/baseball>). Drag coefficient is only a little\n",
    "harder; according to *The Physics of Baseball* (see https://books.google.com/books/about/The_Physics_of_Baseball.html?id=4xE4Ngpk_2EC), the drag coefficient of a baseball is approximately 0.33 (with no units).\n",
    "\n",
    "However, this value *does* depend on velocity. At low velocities it\n",
    "might be as high as 0.5, and at high velocities as low as 0.28.\n",
    "Furthermore, the transition between these values typically happens\n",
    "exactly in the range of velocities we are interested in, between 20 m/s and 40 m/s.\n",
    "\n",
    "Nevertheless, we'll start with a simple model where the drag coefficient does not depend on velocity; as an exercise at the end of the chapter, you can implement a more detailed model and see what effect it has on the results.\n",
    "\n",
    "But first we need a new computational tool, the `Vector` object."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "worthy-wheel",
   "metadata": {},
   "source": [
    "## Vectors\n",
    "\n",
    "Now that we are working in two dimensions, it will be useful to\n",
    "work with *vector quantities*, that is, quantities that represent both a magnitude and a direction. We will use vectors to represent positions, velocities, accelerations, and forces in two and three dimensions.\n",
    "\n",
    "ModSim provides a function called `Vector` that creates a Pandas `Series` that contains the *components* of the vector.\n",
    "In a `Vector` that represents a position in space, the components are the $x$ and $y$ coordinates in 2-D, plus a $z$ coordinate if the `Vector` is in 3-D.\n",
    "\n",
    "You can create a `Vector` by specifying its components. The following\n",
    "`Vector` represents a point 3 units to the right (or east) and 4 units up (or north) from an implicit origin:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "handy-terrain",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = Vector(3, 4)\n",
    "show(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "settled-roommate",
   "metadata": {},
   "source": [
    "You can access the components of a `Vector` by name using the dot\n",
    "operator, like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "vocal-latino",
   "metadata": {},
   "outputs": [],
   "source": [
    "A.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "controversial-shower",
   "metadata": {},
   "outputs": [],
   "source": [
    "A.y"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "earlier-contemporary",
   "metadata": {},
   "source": [
    "You can also access them by index using brackets, like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "digital-channels",
   "metadata": {},
   "outputs": [],
   "source": [
    "A[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "automated-drove",
   "metadata": {},
   "outputs": [],
   "source": [
    "A[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "grave-burst",
   "metadata": {},
   "source": [
    "`Vector` objects support most mathematical operations, including\n",
    "addition:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "conditional-latitude",
   "metadata": {},
   "outputs": [],
   "source": [
    "B = Vector(1, 2)\n",
    "show(A + B)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "charming-reviewer",
   "metadata": {},
   "source": [
    "And subtraction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "encouraging-cabinet",
   "metadata": {},
   "outputs": [],
   "source": [
    "show(A - B)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "combined-command",
   "metadata": {},
   "source": [
    "For the definition and graphical interpretation of these operations, see <http://modsimpy.com/vecops>.\n",
    "\n",
    "We can specify a `Vector` with coordinates `x` and `y`, as in the previous examples.\n",
    "Equivalently, we can specify a `Vector` with a magnitude and angle.\n",
    "\n",
    "*Magnitude* is the length of the vector: if the `Vector` represents a position, magnitude is its distance from the origin; if it represents a velocity, magnitude is its speed.\n",
    "\n",
    "The *angle* of a `Vector` is its direction, expressed as an angle in radians from the positive $x$ axis. In the Cartesian plane, the angle 0 rad is due east, and the angle $\\pi$ rad is due west.\n",
    "\n",
    "ModSim provides functions to compute the magnitude and angle of a `Vector`.  For example, here are the magnitude and angle of `A`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "peripheral-tattoo",
   "metadata": {},
   "outputs": [],
   "source": [
    "mag = vector_mag(A)\n",
    "theta = vector_angle(A)\n",
    "mag, theta"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "great-advice",
   "metadata": {},
   "source": [
    "The magnitude is 5 because the length of `A` is the hypotenuse of a 3-4-5 triangle.\n",
    "\n",
    "The result from `vector_angle` is in radians.\n",
    "Most Python functions, like `sin` and `cos`, work with radians, \n",
    "but many people find it more natural to work with degrees. \n",
    "Fortunately, NumPy provides a function to convert radians to degrees:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "strange-cleaning",
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import rad2deg\n",
    "\n",
    "angle = rad2deg(theta)\n",
    "angle"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "assured-cutting",
   "metadata": {},
   "source": [
    "And a function to convert degrees to radians:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "cellular-community",
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import deg2rad\n",
    "\n",
    "theta = deg2rad(angle)\n",
    "theta"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "responsible-gentleman",
   "metadata": {},
   "source": [
    "To avoid confusion, I'll use the variable name `angle` for a value in degrees and `theta` for a value in radians.\n",
    "\n",
    "If you are given an angle and magnitude, you can make a `Vector` using\n",
    "`pol2cart`, which converts from polar to Cartesian coordinates. For example, here's a new `Vector` with the same angle and magnitude of `A`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "monetary-firmware",
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = pol2cart(theta, mag)\n",
    "C = Vector(x, y)\n",
    "show(C)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "lucky-plastic",
   "metadata": {},
   "source": [
    "Another way to represent the direction of `A` is a *unit vector*,\n",
    "which is a vector with magnitude 1 that points in the same direction as\n",
    "`A`. You can compute a unit vector by dividing a vector by its\n",
    "magnitude:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "explicit-piano",
   "metadata": {},
   "outputs": [],
   "source": [
    "show(A / vector_mag(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "respected-oliver",
   "metadata": {},
   "source": [
    "ModSim provides a function that does the same thing, called `vector_hat` because unit vectors are conventionally decorated with a hat, like this: $\\hat{A}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "relative-republic",
   "metadata": {},
   "outputs": [],
   "source": [
    "A_hat = vector_hat(A)\n",
    "show(A_hat)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "reduced-celebration",
   "metadata": {},
   "source": [
    "Now let's get back to the game."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "similar-local",
   "metadata": {},
   "source": [
    "## Simulating Baseball Flight\n",
    "\n",
    "Let's simulate the flight of a baseball that is batted from home plate\n",
    "at an angle of 45° and initial speed 40 m/s. We'll use the center of home plate as the origin, a horizontal x-axis (parallel to the ground), and a vertical y-axis (perpendicular to the ground). The initial height is 1 m.\n",
    "\n",
    "Here's a `Params` object with the parameters we'll need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "narrative-latest",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "params = Params(\n",
    "    x = 0,          # m\n",
    "    y = 1,          # m\n",
    "    angle = 45,     # degree\n",
    "    speed = 40,  # m / s\n",
    "\n",
    "    mass = 145e-3,    # kg \n",
    "    diameter = 73e-3, # m \n",
    "    C_d = 0.33,       # dimensionless\n",
    "\n",
    "    rho = 1.2,      # kg/m**3\n",
    "    g = 9.8,        # m/s**2\n",
    "    t_end = 10,     # s\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "metric-collins",
   "metadata": {},
   "source": [
    "I got the mass and diameter of the baseball from Wikipedia (see <https://en.wikipedia.org/wiki/Baseball_(ball)>) and the coefficient of drag from *The Physics of Baseball* (see <https://books.google.com/books/about/The_Physics_of_Baseball.html?id=4xE4Ngpk_2EC>):\n",
    "The density of air, `rho`, is based on a temperature of 20 °C at sea level (see <http://modsimpy.com/tempress>). \n",
    "As usual, `g` is acceleration due to gravity.\n",
    "`t_end` is 10 seconds, which is long enough for the ball to land on the ground.\n",
    "\n",
    "The following function uses these quantities to make a `System` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "bored-billy",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from numpy import pi, deg2rad\n",
    "\n",
    "def make_system(params):\n",
    "    \n",
    "    # convert angle to radians\n",
    "    theta = deg2rad(params.angle)\n",
    "    \n",
    "    # compute x and y components of velocity\n",
    "    vx, vy = pol2cart(theta, params.speed)\n",
    "    \n",
    "    # make the initial state\n",
    "    init = State(x=params.x, y=params.y, vx=vx, vy=vy)\n",
    "    \n",
    "    # compute the frontal area\n",
    "    area = pi * (params.diameter/2)**2\n",
    "\n",
    "    return System(params,\n",
    "                  init = init,\n",
    "                  area = area)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "conscious-template",
   "metadata": {},
   "source": [
    "`make_system` uses `deg2rad` to convert `angle` to radians and\n",
    "`pol2cart` to compute the $x$ and $y$ components of the initial\n",
    "velocity.\n",
    "\n",
    "`init` is a `State` object with four state variables:\n",
    "\n",
    "* `x` and `y` are the components of position.\n",
    "\n",
    "* `vx` and `vy` are the components of velocity.\n",
    "\n",
    "When we call `System`, we pass `params` as the first argument, which means that the variables in `params` are copied to the new `System` object.\n",
    "\n",
    "Here's how we make the `System` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "ethical-donna",
   "metadata": {},
   "outputs": [],
   "source": [
    "system = make_system(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "convinced-fellow",
   "metadata": {},
   "source": [
    "And here's the initial `State`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "legitimate-gossip",
   "metadata": {},
   "outputs": [],
   "source": [
    "show(system.init)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "occasional-given",
   "metadata": {},
   "source": [
    "## Drag Force\n",
    "\n",
    "Next we need a function to compute drag force:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "legal-terminal",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def drag_force(V, system):\n",
    "    rho, C_d, area = system.rho, system.C_d, system.area\n",
    "    \n",
    "    mag = rho * vector_mag(V)**2 * C_d * area / 2\n",
    "    direction = -vector_hat(V)\n",
    "    f_drag = mag * direction\n",
    "    return f_drag"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "false-confusion",
   "metadata": {},
   "source": [
    "This function takes `V` as a `Vector` and returns `f_drag` as a\n",
    "`Vector`. \n",
    "\n",
    "* It uses `vector_mag` to compute the magnitude of `V`, \n",
    "and the drag equation to compute the magnitude of the drag force, `mag`.\n",
    "\n",
    "* Then it uses `vector_hat` to compute `direction`, which is a unit vector in the opposite direction of `V`.\n",
    "\n",
    "* Finally, it computes the drag force vector by multiplying `mag` and `direction`.\n",
    "\n",
    "We can test it like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "frank-chick",
   "metadata": {},
   "outputs": [],
   "source": [
    "vx, vy = system.init.vx, system.init.vy\n",
    "V_test = Vector(vx, vy)\n",
    "f_drag = drag_force(V_test, system)\n",
    "show(f_drag)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ultimate-upgrade",
   "metadata": {},
   "source": [
    "The result is a `Vector` that represents the drag force on the baseball, in Newtons, under the initial conditions.\n",
    "\n",
    "Now we can add drag to the slope function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "suitable-salem",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def slope_func(t, state, system):\n",
    "    x, y, vx, vy = state\n",
    "    mass, g = system.mass, system.g\n",
    "    \n",
    "    V = Vector(vx, vy)\n",
    "    a_drag = drag_force(V, system) / mass\n",
    "    a_grav = g * Vector(0, -1)\n",
    "    \n",
    "    A = a_grav + a_drag\n",
    "    \n",
    "    return V.x, V.y, A.x, A.y"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "scheduled-courage",
   "metadata": {},
   "source": [
    "As usual, the parameters of the slope function are a time stamp, a `State` object, and a `System` object. \n",
    "We don't use `t` in this example, but we can't leave it out because when `run_solve_ivp` calls the slope function, it always provides the same arguments, whether they are needed or not.\n",
    "\n",
    "`slope_func` unpacks the `State` object into variables `x`, `y`, `vx`, and `vy`.\n",
    "Then it packs `vx` and `vy` into a `Vector`, which it uses to compute acceleration due to drag, `a_drag`.\n",
    "\n",
    "To represent acceleration due to gravity, it makes a `Vector` with magnitude `g` in the negative $y$ direction.\n",
    "\n",
    "The total acceleration of the baseball, `A`, is the sum of accelerations due to gravity and drag."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "enclosed-favorite",
   "metadata": {},
   "source": [
    "The return value is a sequence that contains:\n",
    "\n",
    "* The components of velocity, `V.x` and `V.y`.\n",
    "\n",
    "* The components of acceleration, `A.x` and `A.y`.\n",
    "\n",
    "These components represent the slope of the state variables, because `V` is the derivative of position and `A` is the derivative of velocity."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "crude-parcel",
   "metadata": {},
   "source": [
    "As always, we can test the slope function by running it with the initial conditions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "closing-simon",
   "metadata": {},
   "outputs": [],
   "source": [
    "slope_func(0, system.init, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adolescent-westminster",
   "metadata": {},
   "source": [
    "Using vectors to represent forces and accelerations makes the code\n",
    "concise, readable, and less error-prone. In particular, when we add\n",
    "`a_grav` and `a_drag`, the directions are likely to be correct, because they are encoded in the `Vector` objects."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "regulated-railway",
   "metadata": {},
   "source": [
    "## Adding an Event Function\n",
    "\n",
    "We're almost ready to run the simulation.  The last thing we need is an event function that stops when the ball hits the ground."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "brief-level",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def event_func(t, state, system):\n",
    "    x, y, vx, vy = state\n",
    "    return y"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "novel-farmer",
   "metadata": {},
   "source": [
    "The event function takes the same parameters as the slope function, and returns the $y$ coordinate of position. When the $y$ coordinate passes through 0, the simulation stops.\n",
    "\n",
    "As we did with `slope_func`, we can test `event_func` with the initial conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "threatened-alberta",
   "metadata": {},
   "outputs": [],
   "source": [
    "event_func(0, system.init, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "unlikely-dressing",
   "metadata": {},
   "source": [
    "Here's how we run the simulation with this event function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "special-background",
   "metadata": {},
   "outputs": [],
   "source": [
    "results, details = run_solve_ivp(system, slope_func,\n",
    "                                 events=event_func)\n",
    "details.message"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "iraqi-appeal",
   "metadata": {},
   "source": [
    "The message indicates that a \"termination event\" occurred; that is, the simulated ball reached the ground.\n",
    "\n",
    "`results` is a `TimeFrame` with one row for each time step and one column for each of the state variables.\n",
    "Here are the last few rows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "prospective-external",
   "metadata": {},
   "outputs": [],
   "source": [
    "results.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dominant-weekly",
   "metadata": {},
   "source": [
    "We can get the flight time like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "medieval-calvin",
   "metadata": {},
   "outputs": [],
   "source": [
    "flight_time = results.index[-1]\n",
    "flight_time"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "convenient-heading",
   "metadata": {},
   "source": [
    "And the final state like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "gorgeous-survey",
   "metadata": {},
   "outputs": [],
   "source": [
    "final_state = results.iloc[-1]\n",
    "show(final_state)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "needed-aruba",
   "metadata": {},
   "source": [
    "The final value of `y` is close to 0, as it should be.  The final value of `x` tells us how far the ball flew, in meters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "failing-bangkok",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_dist = final_state.x\n",
    "x_dist"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "growing-england",
   "metadata": {},
   "source": [
    "We can also get the final velocity, like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "copyrighted-highway",
   "metadata": {},
   "outputs": [],
   "source": [
    "final_V = Vector(final_state.vx, final_state.vy)\n",
    "show(final_V)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "statewide-middle",
   "metadata": {},
   "source": [
    "The magnitude of final velocity is the speed of the ball when it lands."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "structured-adams",
   "metadata": {},
   "outputs": [],
   "source": [
    "vector_mag(final_V)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "vietnamese-diagram",
   "metadata": {},
   "source": [
    "The final speed is about 26 m/s, which is substantially slower than the initial speed, 40 m/s."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "continuous-quick",
   "metadata": {},
   "source": [
    "## Visualizing Trajectories\n",
    "\n",
    "To visualize the results, we can plot the $x$ and $y$ components of position like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "spare-burst",
   "metadata": {},
   "outputs": [],
   "source": [
    "results.x.plot(color='C4')\n",
    "results.y.plot(color='C2', style='--')\n",
    "\n",
    "decorate(xlabel='Time (s)',\n",
    "         ylabel='Position (m)')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "horizontal-bench",
   "metadata": {},
   "source": [
    "As expected, the $x$ component increases as the ball moves away from home plate. The $y$ position climbs initially and then descends, falling to 0 m near 5.0 s.\n",
    "\n",
    "Another way to view the results is to plot the $x$ component on the\n",
    "$x$-axis and the $y$ component on the $y$-axis, so the plotted line follows the trajectory of the ball through the plane:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "dated-browse",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_trajectory(results):\n",
    "    x = results.x\n",
    "    y = results.y\n",
    "    make_series(x, y).plot(label='trajectory')\n",
    "\n",
    "    decorate(xlabel='x position (m)',\n",
    "             ylabel='y position (m)')\n",
    "\n",
    "plot_trajectory(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "certified-synthetic",
   "metadata": {},
   "source": [
    "This way of visualizing the results is called a *trajectory plot* (see <http://modsimpy.com/trajec>).\n",
    "A trajectory plot can be easier to interpret than a time series plot,\n",
    "because it shows what the motion of the projectile would look like (at\n",
    "least from one point of view). Both plots can be useful, but don't get\n",
    "them mixed up! If you are looking at a time series plot and interpreting it as a trajectory, you will be very confused.\n",
    "\n",
    "Notice that the trajectory is not symmetric.\n",
    "With a launch angle of 45°, the landing angle is closer to vertical, about 57° degrees."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "resistant-vegetation",
   "metadata": {},
   "outputs": [],
   "source": [
    "rad2deg(vector_angle(final_V))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cosmetic-aircraft",
   "metadata": {},
   "source": [
    "## Animating the Baseball\n",
    "\n",
    "One of the best ways to visualize the results of a physical model is animation.  If there are problems with the model, animation can make them apparent.\n",
    "\n",
    "The ModSimPy library provides `animate`, which takes as parameters a `TimeSeries` and a draw function.\n",
    "The draw function should take as parameters a time stamp and a `State`.  It should draw a single frame of the animation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "starting-fabric",
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.pyplot import plot\n",
    "\n",
    "xlim = results.x.min(), results.x.max()\n",
    "ylim = results.y.min(), results.y.max()\n",
    "\n",
    "def draw_func(t, state):\n",
    "    plot(state.x, state.y, 'bo')\n",
    "    decorate(xlabel='x position (m)',\n",
    "             ylabel='y position (m)',\n",
    "             xlim=xlim,\n",
    "             ylim=ylim)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sustained-slide",
   "metadata": {},
   "source": [
    "Inside the draw function, you should use `decorate` to set the limits of the $x$ and $y$ axes.\n",
    "Otherwise `matplotlib` auto-scales the axes, which is usually not what you want.\n",
    "\n",
    "Now we can run the animation like this:\n",
    "\n",
    "```\n",
    "animate(results, draw_func)\n",
    "```\n",
    "\n",
    "You can see the results when you run the code from this chapter."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "incomplete-beijing",
   "metadata": {
    "tags": []
   },
   "source": [
    "To run the animation, uncomment the following line of code and run the cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "prescription-boutique",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# animate(results, draw_func)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3bede5d5",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "This chapter introduces `Vector` objects, which we use to represent position, velocity, and acceleration in two dimensions.\n",
    "We also represent forces using vectors, which make it easier to add up forces acting in different directions.\n",
    "\n",
    "Our ODE solver doesn't work with `Vector` objects, so it takes some work to pack and unpack their components.\n",
    "Nevertheless, we were able to run simulations with vectors and display the results.\n",
    "\n",
    "In the next chapter we'll use these simulations to solve an optimization problem."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "lyric-harassment",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "This chapter is available as a Jupyter notebook where you can read the text, run the code, and work on the exercises. \n",
    "You can access the notebooks at <https://allendowney.github.io/ModSimPy/>."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "pressing-retrieval",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    " Run the simulation with and without air resistance.  How wrong would we be if we ignored drag?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "optical-weather",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hint\n",
    "\n",
    "system2 = make_system(params.set(C_d=0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "acknowledged-belgium",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "spatial-ensemble",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "domestic-apparatus",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "correct-pittsburgh",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "generic-shelter",
   "metadata": {},
   "source": [
    "### Exercise 2\n",
    "\n",
    "The baseball stadium in Denver, Colorado is 1,580 meters above sea level, where the density of air is about 1.0 kg / m$^3$.  Compared with the example near sea level, how much farther would a ball travel if hit with the same initial speed and launch angle?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "global-referral",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hint\n",
    "\n",
    "system3 = make_system(params.set(rho=1.0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "appointed-sugar",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "specialized-mediterranean",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "shaped-paragraph",
   "metadata": {},
   "source": [
    "### Exercise 3\n",
    "\n",
    " The model so far is based on the assumption that coefficient of drag does not depend on velocity, but in reality it does.  The following figure, from Adair, *The Physics of Baseball*, shows coefficient of drag as a function of velocity (see <https://books.google.com/books/about/The_Physics_of_Baseball.html?id=4xE4Ngpk_2EC>).\n",
    "\n",
    "![Graph of drag coefficient versus velocity](https://github.com/AllenDowney/ModSimPy/raw/master/figs/baseball_drag.png)\n",
    "\n",
    "I used an online graph digitizer (<https://automeris.io/WebPlotDigitizer>) to extract the data and save it in a CSV file.  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "german-niagara",
   "metadata": {
    "tags": []
   },
   "source": [
    "The following cell downloads the data file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "directed-moisture",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "download('https://github.com/AllenDowney/ModSim/raw/main/data/' +\n",
    "         'baseball_drag.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "imported-chosen",
   "metadata": {
    "tags": []
   },
   "source": [
    "We can use Pandas to read it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "fuzzy-register",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from pandas import read_csv\n",
    "\n",
    "baseball_drag = read_csv('baseball_drag.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "progressive-buyer",
   "metadata": {
    "tags": []
   },
   "source": [
    "Here are the first few rows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "returning-fellowship",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "baseball_drag.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bronze-acoustic",
   "metadata": {
    "tags": []
   },
   "source": [
    "I'll use Pint to convert miles per hour to meters per second."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "reasonable-swaziland",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "mph_to_mps = (1 * units.mph).to(units.m/units.s).magnitude\n",
    "speed = baseball_drag['Velocity in mph'] * mph_to_mps"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "industrial-architecture",
   "metadata": {
    "tags": []
   },
   "source": [
    "I'll put the results in a `Series`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "attached-shower",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "C_d_series = make_series(speed, baseball_drag['Drag coefficient'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "amazing-horse",
   "metadata": {},
   "source": [
    "Here's what it looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "id": "christian-camcorder",
   "metadata": {},
   "outputs": [],
   "source": [
    "C_d_series.plot(label='$C_d$')\n",
    "decorate(xlabel='Speed (m/s)', \n",
    "         ylabel='Coefficient of drag')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "computational-alloy",
   "metadata": {
    "tags": []
   },
   "source": [
    "And, for use in the slope function, we can make a function that interpolates the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "heated-belfast",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "drag_interp = interpolate(C_d_series)\n",
    "drag_interp(30)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "reverse-shock",
   "metadata": {},
   "source": [
    "Modify the model to include the dependence of `C_d` on velocity, and see how much it affects the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "engaged-provision",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "directed-fiber",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "framed-dealer",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "id": "accomplished-elizabeth",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "going-techno",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "brief-saying",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "spare-pregnancy",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "id": "catholic-staff",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "id": "broad-sequence",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "strong-design",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
